BOUNDARY DRIVEN WAVEGUIDE ARRAYS: 
SUPRATRANSMISSION AND SADDLE-NODE BIFURCATION* 

HADI SUSANTOt 

Abstract. In this report, we consider a semi-infinite discrete nonlinear Schrodinger equation 
driven at one edge by a driving force. The equation models the dynamics of coupled waveguide 
arrays. When the frequency of the forcing is in the allowed-band of the system, there will be a linear 
transmission of energy through the lattice. Yet, if the frequency is in the upper forbidden band, then 
there is a critical driving amplitude for a nonlinear tunneling, which is called supratransmission, 
of energy to occur. Here, we discuss mathematically the mechanism and the source of the supra- 
transmission. By analyzing the existence and the stability of the rapidly decaying static discrete 
solitons of the system, we show rigorously that two of the static solitons emerge and disappear in a 
saddle-node bifurcation at a critical driving amplitude. One of the emerging solitons is always stable 
in its existence region and the other is always unstable. We argue that the critical amplitude for 
supratransmission is then the same as the critical driving amplitude of the saddle-node bifurcation. 
We consider as well the case of the forcing frequency in the lower forbidden band. It is discussed 
briefly that there is no supratransmission because in this case there is only one rapidly decaying 
static soliton that exists and is stable for any driving amplitude. 
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1. Introduction. An exotic nonlinear phenomenon has been discovered recently 
by Geniet and Leon [1] in a semi-infinite chain of coupled oscillators driven at one 
edge by a time periodic forcing. Energy excitations will propagate through the chain 
if the driving frequency is in the allowed-band of the discrete system. It is natural be- 
cause of the system's dispersion relation. In contrast, it would be expected that if the 
forcing frequency is in the band-gap, then there would be no energy flow. Yet, Geniet 
and Leon [1] show theoretically and experimentally that there is a definite driving 
amplitude threshold above which a sudden energy flow takes place. This phenomenon 
is called nonlinear supratransmission T . An exciting independent work on a modified 
Klein-Gordon equation describing the Josephson phase of layered high-Tc supercon- 
ductors shows the presence of the same phenomenon [3]. Promising technological 
applications employing supratransmission have been proposed as well accompanying 
these findings, such as binary signal transmissions of information [4] and terahertz 
frequency selection devices [5]. 

In [6] Khomeriki considers boundary driven coupled optical focusing waveguide 
arrays described by 

(LI) = -^n+l - i'n-l - ll^nl^^n, ^0 = Ae'^' , 

OZ 

with 7 > and n = 1,2,.... Here, Tpn is the electromagnetic wave amplitude in 
the nth guiding core, z is the propagation variable, A is the propagation constant 
or the driving frequency, and 7 represents the nonlinearity coefficient which is taken 
to be 7 = 2 in this report. This model can also be considered as a slow modulation 
wave approximation to the discrete sine-Gordon equation [7j. Similar to the nonlinear 
band-gap tunneling observed by Geniet and Leon [T], it is reported that there is a 
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Fig . 1.1. Three dimensional plots of time evolution of the boundary driven waveguide arrays 
Eq. ifj.iD . When the driving frequency —2 < A < 2 is in the allowed band, any driving amplitude 
will lead to an energy flow to remote sites (top left panel). If A is in the upper forbidden band and 
A is small enough, the boundary will excite a couple of arrays only (top right panel). Yet, there 
is a critical threshold amplitude ylth(A) above which there is nonlinear forbidden band tunneling 
indicated by the released of a train of discrete solitons (bottom left panel). A quantitatively different 
behavior of supratransmission occurs when the driving frequency is large enough as is shown in the 
bottom right panel. 

critical threshold ylth(A) for supratransmission when the propagation constant A is 
in the forbidden band A > 2 [6J . 

In Fig. II. 1[ we present numerical simulations of the dynamics of Eq. Fol- 
lowing [B] , the driving is turned on adiabatically to avoid the appearance of an initial 
shock by assuming the form 

A = A{1- exp(-z/T)), 

where we omit the breve henceforth. In the following figures, we take t — 50 and 
apply a linearly increasing damping to the last 20 sites to suppress edge reflection. 

Presented in the top left panel of Fig. 11.11 is a three dimensional plot of time 
evolution of Eq. (|l.ip when the driving frequency is in the allowed band —2 < A < 2. 
A small driving amplitude will excite all the sites. On the other hand, if the driving 
frequency is in the upper band A > 2, a small A will only excite several neighboring 
sites as is shown in the top right panel of Fig. II. II Yet, if the driving amplitude is large 
enough, then a train of 'traveling' discrete solitons can be released [6] (see bottom left 
panel of the same figure). This flow of energy is the so-called supratransmission or 
nonlinear forbidden band tunneling and we call the minimum A for supratransmission 
to occur a critical threshold Ath- The word traveling is in quotes because Eq. (|l.ip 
does not admit a genuine one (see [9] and references therein). If one waits long 
enough, the gap-solitons will be trapped by the lattice. Khomeriki [6] also notices 
an immediate trapping when the driving frequency A is relatively large as is shown 
in the bottom right panel of Fig. 11.11 In this regime, the corresponding discrete gap 
solitons are highly localized. 

An analytical approximation of Ath (A) hi the limit 0<A — 2<Clis given by [B] 

(1.2) Ath(A) = VA^. 

Remark 1.1. Equation U.l\) is symmetric with respect to the transformation 
A — !■ ~A and ipn ~> This means that there is also a critical amplitude —Ath (A) 

if one applies A < such that for A < —Ath (A) < 0, a nonlinear forbidden band 
tunneling will occur. Equation il.l]) is also symmetric with respect to the transfor- 
mation A — ^ —A, — > (~l)"V'n; o^*^ 7 ^ ^7- Therefore, the same phenomenon 
can be observed in defocusing waveguide arrays 7 < 0. Due to the transformation, the 
only difference of defocusing arrays from the self-focusing ones is that there will be a 
TT phase- difference between neighboring lattices. 

It is presented in [6j that the numerical result for the threshold amplitude deviates 
rapidly from the approximation (|1.2p . It is because Eq. (|1.2[) is actually the amplitude- 
temporal frequency relation of the continuous nonlinear Schrodinger (NLS) equation's 
solitons. The relation has been phase-shifted properly due to some transformation, i.e. 
tpn ipn exp{2iz). Applying the transformation to (|l.ip will take it to a normalized 
standard finite difference approximation of the continuous NLS equation. 
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Remembering the aforementioned promising applications of nonlinear tunneling, 
it is therefore of interest to obtain an approximation of the threshold amplitude in 
the other limit A — 2 ^ 1. It is one of the aims of the present report. The other 
aim is to understand mathematically the mechanism of the nonlinear tunneling. It is 
mentioned, but not rigorously proved, ^ that the supratransmission happens because 
of the emergence of two solutions at the critical driving amplitude, i.e. a saddle- 
node bifurcation. If it is the case, this then means that the threshold amplitude 
is not necessarily the amplitude of the corresponding fundamental soliton Eq. (|1.2p . 
Understanding the source of supratransmission also will allow us to explain, e.g., the 
reason why there is no threshold amplitude for nonlinear tunneling when the driving 
frequency is in the lower forbidden band A < — 2. 

Nonetheless, one may question the relevance of our first aim, as supratransmission 
is quickly trapped by the lattices for large A. Even though our analysis may be not 
immediately applicable to the present case, the aim is still of relevance. There are sev- 
eral experimentally realizable discrete equations that support 'traveling' solitons in a 
parameter region where the gap solitons arc highly localized. As a particular example 
is the discrete Schrodinger equation with saturable nonlinearity in the large nonlinear- 
ity coefficient regime flOJ. We have observed supratransmission in this equation and 
have successfully applied our analysis presented herein to obtain an approximation 
to the threshold amplitude [TT]. Later on in this paper, we also conjecture that our 
analytically obtained approximation, presented in terms of a power series expansion, 
may well be convergent uniformly in the region of interest, i.e., A > 2. Moreover, 
the mathematical procedure presented herein can also be applied as an alternative 
method to analyze the bistability effect considered, e.g., in [12]. We might even con- 
sider it simpler and more appropriate as the analysis can then be done solely in its 
discrete set-up, with no necessity of approximating the problem with its continuous 
counterpart [12]. 

In this study, we will show that the supratransmission is indeed related to saddle- 
node bifurcations. To mathematically prove this, our strategy is as follows. We will 
first prove the existence of a mode bifurcating from the constant solution ■(/'„ = 
due to the driving site. We will also show that there is a singular mode bifurcating 
from infinity. We will then demonstrate that these two modes collide in a saddle-node 
bifurcation by developing an asymptotic analysis in the range of A large. Such an 
analysis is doable in that regime because the modes are highly localized. The final 
step to show that the critical amplitude is the same as the threshold amplitude for 
supratransmission is to prove that the mode bifurcating from zero state is stable, all 
the way on its existence region. Using this result, then we can derive an approximation 
of the threshold amplitude in terms of a power series expansion that can be calculated 
to any order. Numerical computations will be presented as well to compare our 
analytical results. 

Our paper will be outlined in the following structures. In Subsection 2.1, we 
present our asymptotic analysis for the existence of monotonically decaying static 
solutions of Eq. (|l.ip . The next subsection will contain our study on the stability 
analysis of solutions discussed in the preceding subsection. Using the same procedures, 
we then briefly discuss in Subsection 2.3 that there is no supratransmission in the 
case of A < —2 as there is no bifurcation occurring in this regime. Then, we compare 
our analytical findings with the results of numerical computations in Subsection 2.4. 
Finally, we summarize our findings and present our conclusions in the last section. 

2. Existence and stability analysis of rapidly decaying discrete solitons. 
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2.1. Existence analysis. Stationary solutions of Eq. (jl.ip are sought in the 
form of V'nl-z) = 0ne*^^j where 0„ is a real valued function. This ansatz can be 
applied as one would naturally expect that all the sites will be excited with the same 
frequency as the driving frequency. Since we are interested in the large propagation 
constant A, we scale A — > 1 and define e = 1/A. Hence, we consider |e| ^ 1. Static 
equation of Eq. is then given by 



with (jjQ = A. 

When |e| is small enough, apart from the boundary, the leading order solution of 
(pn would formally satisfy 



arrays are almost uncoupled and indicates that solutions of Eq. (|2.ip can be expressed 
in terms of an asymptotic or a perturbation expansion in e. It also says that when 
we consider finitely long waveguide arrays, i.e. n = 1,2, . . . , N , Eq. (|2.ip can have 
at most 3^ solutions. Yet, only some of them are related to the nonlinear tunneling 
phenomenon presented in Fig. 11.11 We are especially interested in solutions with a 
magnitude that is monotonically decaying with the property |(/)„| ^ as n — > oo. This 
consideration is based on the fact that when the driving frequency is in the forbidden 
band and the driving amplitude is below the critical threshold, the solution profile is 
monotonically decaying as n ^ oo (see the top right panel of Fig. II. ip . Moreover, we 
only need to consider particularly a family of rapidly decaying discrete solitons which 
is defined as follows. 

Definition 2.1. Let cpn — ^2^=0 o,n,k'&k{() be a solution of \2.1\) . where ri e Z+, 
'dk{e) is an asymptotic sequence and i)k{e) — o(i?fc_i(e)) for e — > 0. 0„ is a rapidly 
decaying discrete soliton of Eq. 112. 1\) if\(f>n\ is a monotonically decreasing function to 
as n oo with a property that to the leading order O{'do) only the first lattice site 
is non-zero, i.e. ai^ ^ and a„^o = 0, n 7^ 1. 

As an example of this definition, let us consider the following solution of (|2.ip 



This solution is obtained from the expansion: (\)\ = — l/Y^e7 + ai_i\/e + ai.2e + . . ■, 
02 = l/\/e7 + a-i.wf^ + 02,26 + . . ., (/>3 = + 03^1 + • • •, and 0„ = + . . . for 
n > 3. Substituting the ansatz to Eq. (|2.ip will yield polynomials in e. Equating the 
coefficients of the polynomials for all orders of e to zero will yield equations for ak,i 
that have to be solved simultaneously to obtain Eq. p.3p . 

It is clear that the profile of |$o("-7^)| (|2.3p is monotonically decaying in n. 
However, this solution is not rapidly decaying as to the leading order, i.e. 0(l/^/e), 



The existence of rapidly decaying solutions of ()2.ip when A = 0(1) is guaranteed 
by the following theorem. 

Theorem 2.2. Let A he of 0{l). Then for e positive and small there are three 
rapidly decaying discrete solitons of the static equation Denoted by j = 



(2.1) 





(2.3) 




|$o(2,^)| = |$o(l,^)|^0. 



BOUNDARY DRIVEN WAVEGUIDE ARRAYS 



5 



1,2,3, the solitons are given by 



(2.4) 



$i(n, A) = 



(2.5) 



(2.6) 



$3(n,A) = 



A-- — 
2 2V7 

0(e2), n = 2 
-0(6^), n = : 



+ 0(e2), n = 1 



V7 

£3/2 

+ 0(6 ), otherwise, 

Ae + Ae^ + 0{e'), n. = l 
Ae^ +Oie^),n^2 
Ae3 + 0(e5), n = 3 
+ C'(e''), otherwise, 

A- + -^ + C'(e2),n=l 



67 

g3/2 



2V7 



+ 0(e2), n 



— +0(e^), n = 3 
V7 

+ 0(6 ), otherwise. 



Proof. Because we are looking for rapidly decaying solitons, to the leading order 
Eq. ()2.1|) can be represented by 



(2.7) 



bi+eA + 7601 = 0. 



Equation (|2.7p is a cubic equation similar to Eq. (|2.2p . also with three roots. However, 
as e — > 0, (|2.7p reduces to a linear equation 01 = with only a single root. Therefore, 
finding the roots of the equation is a singular perturbation problem. Following, e.g., 
[14j (see Example 3 of Section 2.1 and 2.2), one will obtain the roots of (|2.7p . i.e. (j>i = 
Ae + . . . and 0i = ±l/y/je+ .... This concludes that there are three rapidly decaying 
solutions of (|2.ip . In the following, let us name the solitons ^j{n, A), j = 1, 2, 3, with 
$i(l,v4) = 1/Vi7+ *2(l,v4) = e^+ and $3(1,^) = -1/^67 + The 
existence of ^j{n,A) for Eq. (|2.ip follows immediately from the Implicit Function 
Theorem (see, e.g., [I3j) since F is differentiable and the Jacobian matrix of problem 
(|2.ip DF{(t),0) is invertible. Explicit calculations to obtain (|2.4p - (|2.6p can be done 
similarly following the derivation of (|2.3p . □ 

If one compares the above theorem and the top right panel of Fig. 11.11 it can be 
recognized immediately that the solution observed in the panel in the limit 2; ^ 00 is 
nothing else but |$2(»^,^)|- 

One still can obtain the existence of the above rapidly decaying solutions even 
when A ^ 1 as stated in the following theorem. 

Theorem 2.3. Let A be scaled to A = i/e^/^ A < 2/^/2Tj. 



$0 

3 



A 



(2.8) 



$,(n,A) 



37$f 



1 



-0(e3/2), 



$°V^ + 0(e3/2), „ = 2 
+ 0(e^/^), otherwise, 
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with <f>^ given by 



(2.9) $^ = < 



2 



2 



cos 



47r 1 

— + - arccos 
2n 1 

— + - arccos 



-Ay/TFf 



J =3. 



Moreover, if we write A = 2/^(277e^) - with A > l/y^, then ^j, j = 1,2, 

can be written as 



(2.10) 



1 



87 87 



■ + 0(6^/2), n = l 



3j 



0{e- 



3/2N, 



n = 2 



0{e^^^), otherwise. 



Proof. As we are interested in the case of A 1, we first scale A = Aje"^"^ and 
correspondingly write $j(n, A) — ^°{n, A)/y/e+ . . j ~ 1, 2, 3, with $°(n, A) — for 
all 1 < ?i e Z+. Substituting the expansion to Eq. (|2.ip and identifying coefficients 
for power series of 0{l/y/e) yields the following cubic equation for <&°(1, A) — i.e. 



(2.11) 



G($o) :=-$° + i + 7(<i>°)-^=0. 



Equation (|2.1ip cannot be solved perturbatively to obtain the roots $^ as before 
as all the terms in (|2.1ip are of the same order. Therefore, we need the following 
lemma on cubic equations. 

Lemma 2.4. Consider the following polynomial equation 



g{x) = ax^ + bx^ + cx + a, b,c,d d 



Let 



Y = g{X), 

da 



2aV^ 



Sac 



9a2 



-arccos( — ). 



IfY^ < , then the cubic equation has three distinct real roots given by 



(2.12) 
(2.18) 
(2.14) 



xi — X + 2vcos9, 
X2^X + 2ucos(47r/8 + 0), 
X3. = X + 2i;cos(27r/8 + 6'), 



where 



Xi> X2> x^. 
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When — h? , two of the roots which are neighboring to each other, i.e. Xi and 
or Xi and x^, will collide in a saddle-node bifurcation and disappear when > , 
i.e. the cubic equation then only has a single real root. 

Proof. See [T5|. The expression of the cubic polynomial roots (|2.12l) - (|2.14p is 
using Nickalls' geometric representation [16]. □ 

According to Lemma [2^ Eq. (|2.1ip has geometric representation parameters: 

1 2 1 —Y 

X = Q, Y = A, V ~ , h — , 6 = — arccosf ), 

^/3^ V277 3 ^ h 



from which we can conclude that (|2.1ip has three real roots when A < 2/y/27'y. 
The roots of ^J^, i.e. Eqs. (EH), arc obtained using Eqs. ((TT^ - fTTO)) . Then, the 
continuation of $j can be obtained immediately using the Implicit Function Theorem. 

It is then straightforward to calculate that when A = 2/y/TFf, 2 — 
2cos0 = 2cos(47r/3 + 0) — 1. Hence, we know that ^i{n,A) collides with ^2{n,A) 
in a saddle-node bifurcation. 

For the value of A close to the occurrence of the saddlc-nodc bifurcation, we 
write A = 2/^/27je^ — A^/e. In this case, the Implicit Function Theorem cannot be 
immediately employed to prove the existence of $1 and $2 as we need a bound for A. 

First, we substitute to the steady state equation (|2.ip <i>j = <i>°(n,A)/y^ + 
\/e$](n,A), j = 1,2, with ^°{n,A) = for n = 1 and otherwise. This 

then gives the following equations: 

Gi ($] , e) := $] (2, A) - 1 + ($] (1, ^)) ' + £7 ($1 (1, A)) ' = 0, 
G2($],e) := -$](2,A) + e$](3, A) + -|= + e<i>](l, A) + e^^ (<i>](l, A))' = 0, 
G„($], e) := A) + e +l,A) + - 1, A) + e-/<^>]{n, A)^) =0, n ^ 1, 2. 

Taking e = 0, the above equations give us 

$](1,A) = ±^' ^ 

^]i2,A) = 
*](n,A)=0, n^l,2. 

Note that the zb-solutions collide for A = 1/^/37. Because the linearization I?G($j, 0) 

is invertible for A > 1/ -y/37, the Implicit Function Theorem can be applied again and 
we have the existence of rapidly decaying solitons = ^^{n, A)/ ^/e + y/e^^{n, A), 
j = l,2. □ 

With this theorem, we then have shown that $1 collides in a saddle- node bifur- 
cation with $2- Yet, we cannot directly claim that this is the source of the supra- 
transmission observed in Fig. Il.ll before we show and discuss the stability of the two 
solitons. 

2.2. Stability analysis. After discussing the existence of rapidly decaying soli- 
tons of Eq. ()2.ip . next we study their stability. If (/>„, n — 1,2,..., is a solution 
of (|2.ip . then the linear spectral stability of 0„ can be obtained by substituting the 
ansatz Vn {(f>n + (5Ke*^^ -f- uJ;:e-*^^])e'^^ with A G C, (u„, w„) £ C^, and n e Z+ 
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into Eq. (|l.ip . Linearizing the equation to 0{S), we obtain the foUowing eigenvalue 
problem 

(2.15) Xe f ) - f )+^( "'')+^^( 
with 

vo \ f 0\ f 1 

wo)^[o)^ ^=(,0 -1 

V -n(l>n 1 -2e7|0„|2 y 

where we have scaled A ^ 1. 

The natural domain for C = {ea C ea) is L^(C). We call A an eigenvalue of 
C if there is a function {w„}„gz+; {w,i}n<£Z+ G L^(C) which satisfies (|2.15p . Since C 
depends smoothly on A, the eigenvalues of £ will depend smoothly on A, too. 0„ is 
linearly stable if the imaginary part of A is zero, i.e. Im(A) = 0. 

The continuous spectrum is obtained by substituting 

to Eq. (|2.15p from which we will obtain 

eA = ±2e cos fc ^ 1. 

Thus, the continuous spectrum of solutions under investigation is the range 

(2.16) Ae ^-i-2,-i + 2^ andAe Q - 2, ^ + 2 

As the continuous spectrum lies in the real axis, the stability of the solutions is only 
determined by the discrete spectrum, i.e. eigenvalues. For the solutions given in 
Theorems 12.21 and 12. 3[ we have the following stability results. 

Theorem 2.5. For small driving amplitude A — 0(1), the various rapidly de- 
caying discrete solitons have the following properties: 

1. the discrete soliton <i>i is unstable. It has a single imaginary eigenvalue. 

2. the soliton $2 is strictly stable as the soliton has no discrete eigenvalues. 

3. the discrete soliton $3 is stable. It has a single real eigenvalue. 

Proof. We are looking for eigenvectors that are also rapidly decaying. Therefore, 
the eigenvalue problem Eq. (|2.15p to the leading order can be approximated by the 
linear eigenvalue problem 

Aef "1 \^C' "1 



Wi I \ Wi 

which gives the following approximate eigenvalues 

1 



(2.17) A = ±-V3(e7(/.„2) -Ae-fc^J + l. 

In the above expression, we have taken into account the fact that 0„ e 
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For the stability of <i>i and $3, substitute 0i = «f>j(l, A), j = 1, 3 into Eq. (|2.17p . 
Taking the series expansion of the expression gives the foUowing eigenvalue A for 
$j(n, A), i.e. 



(2.18) A = 



Because the eigenvalue of ^i{n,A) has a non-zero imaginary part, we conclude that 
to the leading order $1 is unstable, as opposed to $3. 

As for — <i>2(l,v4), the series expansion of Eq. (|2.17p gives 

(2.19) A = l/e + 0(e2). 

Because A is inside the continuous spectrum (j2.16p . then our assumption that the 
eigenfunction is rapidly decaying is not justified. Nonetheless, we know that $2 bifur- 
cates from a uniform solution 0„ = which is stable. Because L depends smoothly 
on A, we then can conclude that <i>i has no eigenvalue. □ 

When the driving amplitude is large, we also have the following theorem 
Theorem 2.6. For large driving amplitude A — Ae~'^l'^ , the various rapidly 
decaying discrete solitons have the following properties: 

1. the discrete soliton $1 is unstable with a single imaginary eigenvalue. 

2. the soliton $2 strictly stable with a single real eigenvalue. 

3. the discrete soliton $3 is in general stable with a single eigenvalue, except in 
a finite interval where our asymptotic analysis is inconclusive. 

To the leading order, the eigenvalue of the three solitons is given by 



(2.20) \ = K/e+ -. ^ ^ ^ ^- + K\+0{^fe) 




with K = \I3 (7*°^) - 47$°^ -I- 1. Moreover, by writing A = 2/ y/27je^ - A^/e, the 
eigenvalue o/$i,2 is given by 

2 



(2.21) A = 



31/4 Vi\ 




with the 'minus' sign for the eigenvalue o/$i and the 'plus' sign for $2- 

Proof. The proof of Theorem 12.61 is similar to the proof of Theorem 12.51 The 
stability result of <I>3 cannot be deduced immediately because the expression of $3 
is not trivial. The presence of a finite interval where our asymptotic analysis is 
inconclusive cannot be seen clearly. It is inconclusive because there is a range of A in 
which A is in the domain of the continuous spectrum (|2.16p . A numerical proof will 
be presented in the following section. □ 

2.3. Analysis for the case of A < —2. We omit the details and the rigorous 
proof, but it can be shown that for A < — 2, there is only one rapidly decaying soliton 
which is stable for any driving amplitude. The idea is as follows. 

Instead of Eq. (|2.ip , consider 



(2.22) 



0n = -e((?!>n+l + - 7e0n^: 4'0 = A, 
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where we again have scaled > A — > 1 and define e — 1 / 1 A | . For a rapidly decaying 
solution, the leading order equation of (I2.22p is then given by 



It is clear that / ±oo as </> — > ±oo. Yet, / has no critical point, i.e. df/dcfii > 
0. Therefore, one can conclude that / is a monotonically increasing function which 
intersects the horizontal axis once, i.e. / has one real root. The stability of this rapidly 
decaying solution might be determined immediately following Theorem 12.51 and 12.61 
Our numerics, which are not presented here, show that when A = 0(1) the solution 
is stable with no discrete spectrum and when A = 0{l/€^^^) there is an eigenvalue 
bifurcating from the upper edge of the continuous spectrum. Hence, the soliton is 
stable all the way to ^ ^ oo which explains why there is no supratransmission for 



2.4. Numerical results. To accompany our analytical results, we have used 
numerical calculations. For that purpose, we have made a continuation program 
based on a Newton iteration technique to obtain stationary rapidly decaying discrete 
solitons of Eq. (12. 1|) and an eigenvalue problem solver to solve ()2.15|1 in MATLAB. 
Throughout the subsection, we consider in particular A = 10. Even though there is 
no prominent supratransmission of energy for this value of A, it is taken solely as 
an example to show that especially in the regime of A large, our asymptotic analysis 
explains the problem well. It will be shown below that, e.g., even using the first 
two terms of the approximate threshold amplitude, our analytical result is already 
relatively in agreement with the numerical results. 

We summarize our results and discussions for the existence and the stability of 
<I>i and $2 in Fig- 12.11 At the top left panel of the figure, we present the existence of 
<I>j, j — 1, 2, represented by the solution of the first site, where the upper and lower 
branch corresponds to the existence curve of $i and $2, respectively. Presented in 
solid line is the numerical results. Our analytical result Eqs. (|2.4p and (|2.5p which 
is supposed to be vahd when A = 0(1) is depicted as dash-dotted line. As for the 
analytical approximations for A = 0(l/e^/^), i.e. Eqs. ()2.8p and (I2.10p . they are 
presented as dotted and dashed line, respectively. It is interesting to note that Fig. 
l2.1l shows clearly a good agreement between our analytical and the numerical results. 

Top right panel of Fig. 12.11 presents the comparison between the critical ampli- 
tude ^th(A) calculated numerically from Eq. (j2.ip and our approximation Ath(A) = 
2/ ^/2T-fe^ — (see Theorem 12. 3p which are presented in solid and dashed 

line, respectively. The numerical results were also checked against the full dynamics 
of the original problem Eq. (jl.ip . where an agreement is obtained as it should be 
provided that r is large enough. Note the good agreement when A ^ 1. As a com- 
parison with the analytical approximation obtained by Khomeriki [6] , we also present 
Ath(A) = y/ A — 2 in dash-dotted line. 

It is interesting to note that in the limit A ^ 2 our analytical approximation does 
not diverge. As is shown in the inset of the top right panel figure, the difference of 
the approximate value of the threshold amplitude and the numerical result at A = 2 
is about 50%. Using the same method presented in the preceding sections, we obtain 
that the first three terms of the approximation of Ath(A) are actually given by 



(2.23) 




A < -2. 



(2.24) 
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Fig. 2.1. Presented is the comparison between the numerically obtained results and the analyti- 
cal calculations presented in Section 2. Top left panel: the existence curve of <I>i and <I>2 represented 
by the solution of the first site, where the upper and lower branch corresponds to the existence curve 
of^l and <I>2, respectively. Top right panel: the threshold amplitude Atjj as a function of the propa- 
gation constant A. Bottom panels: the critical eigenvalue of (left) and $2 (right) as a function 
of the driving amplitude A. Shaded region in the bottom right panel shows the region for the con- 
tinuous spectrum. Analytical approximations calculated in Section 2 are also presented as dashed, 
dotted, and dash-dotted lines (see the text). 

The plot of this curve is depicted in the same panel as dotted line where one can 
see that the difference now has decreased by about 10%. This then motivates us to 
question whether the infinite power series of the approximate threshold amplitude 
Ath(A) is actually convergent uniformly to the critical amplitude curve. Considering 
the fact that the region of interest is on < e < 1/2 and the coefficients of the 
power series are so far bounded, the answer might well be affirmative. Yet, this 
question is out of the scope of the present report and it is therefore addressed for 
future investigations. 

After presenting the numerical and the analytical results for the existence of $1 
and $3, next we consider the stability of the solitons. Bottom panels of Fig. 12.11 
present the comparison between the results. The left panel shows imaginary part of 
the critical eigenvalue of $1 as a function of A in its existence region. It is clear 
that the soliton is always unstable. The right panel presents the eigenvalue of ^2 as a 
function of the driving amplitude where one can see that the soliton is always stable, as 
opposed to $1. Our analytical approximations (|2.18[) . (|2.20p . and (|2.21D are presented 
as well in the two panels as dash-dotted, dotted, and dashed line, respectively. It is 
also interesting to note that as is predicted by Theorem l2.21 $2 has no eigenvalue when 
A is small. Our analytical approximation (I2.2ip predicts very well the appearance of 
the eigenvalue of $2- 

Because it is known that $1 is unstable in its entire existence region, it is of 
interests to see how the dynamics concerning the instability. In Fig. (j2.2p we present 
the evolution of $1 for a parameter value A = 8.46 {A is already at this value from the 
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Fig. 2.2. The instability dynamics o/ $1 for A = 8.46 and A = 10. It is presented that even 
with a tiny different perturbation the dynamics can be significantly different (see the text). 

beginning z = 0, as opposed to Fig. 11.11 where ^ is in the beginning and gradually 
increases to a constant). The top left panel presents the dynamics of $1 with the 
initial condition = 0) = ^i{n,A) — 10""*. The initial condition ^i{n,A) is 

obtained numerically from Eq. I|2.ip . The top right panel depicts the behavior of 
the first site in time where one can see that the instability manifests in the form of 
soliton's oscillations. Interestingly, if we start with an initial condition of the form 
ipn{z = 0) = $i(n, A) + 10~^, the solution has a similar instability behavior but with 
a different oscillation maximum. The dynamics is presented in the bottom panels of 
Fig. 12.21 It is important to note that with such a small change, the dynamics can be 
significantly different. This duality therefore might be employed as a small intensity 
light detector similar to the proposal of [T^ . 

We have analyzed as well numerically the existence and the stability of the soliton 
$3. We summarize our results in Fig. 12.31 The numerical result for the existence of 
the soliton is shown in the top left panel of the figure. Our analytical approximations 
(|2.6p and (|2.8p are shown in dash-dotted and dotted lines, respectively, where one can 
see the good agreement between the numerical and the analytical result. 

After studying the existence of the discrete soliton, we next present our stability 
analysis of the soliton. Shown in the bottom panels of Fig. 12.31 is the numerically 
obtained critical eigenvalue of $3 as a function of A. In the bottom left panel is the real 
part of the critical eigenvalue. It is clear that when ^ = 0, the eigenvalue is a double 
eigenvalue at zero. As soon as A is increased, the zero eigenvalue bifurcates along 
the real line. At a critical driving amplitude, the eigenvalue collides with the lower 
boundary of the continuous spectrum. The result of the collision is the bifurcation 
of the eigenvalue into the complex plane resulting in an eigenvalue with nonzero 
imaginary part. In the bottom right panel, we present the trajectory of the eigenvalue 
in the complex plane as A is increased. One can then see that there is also another 
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Fig. 2.3. Similar to Fig. \2.l\ but for the soliton $3. Top left panel shows the numerical results 
for the existence of ^3 vs. the driving amplitude A (solid line). Presented is the value of the solution 
at the first site, i.e. <&3(1, A). Bottom left panel presents the stability of the soliton. The red solid 
line that separates the black solid line indicates that the soliton is unstable in this region. The 
behavior of the critical eigenvalue in the complex plane is depicted in the bottom right panel. In the 
panel, the parametric variable is the driving parameter A. The top right panel shows the dynamics 
of the soliton when it is unstable. See the text for the analytical approximation curves. 

critical amplitude above which the eigenvalue becomes real again, i.e. the soliton 
becomes stable. In the region where the imaginary part is nonzero, we depict the 
curve in the bottom left panel of Fig. 12.31 in solid red line. We also compare it 
with our analytical approximations Eqs. (|2.18p and (|2.20|) . that are shown in dash- 
dotted and dotted line, respectively. In Theorem 12.61 it is stated that our analytical 
approximation is inconclusive for the case of A large. As can be seen from Fig. 12.31 
our analytical approximation Eq. (j2.20p is always real. It is because when the real 
part of the eigenvalue is in the region of the continuous spectrum, our assumption 
that the eigenfunction is fastly decreasing is not justified. 

It is then interesting to see the dynamics of the instability. In the top right panel, 
we depict the evolution of an unstable discrete soliton of type $3. The parameter 
values are depicted in the figure. The setup for the driving amplitude is similar to 
the setup of Fig. O 

Regarding the involvement of $3 in the dynamics of the driven boundary waveg- 
uides (jl.ip (see Fig. II. ip . it is not clear whether when $2 disappears it evolves into 

$3- 

3. Conclusions. We have analyzed mathematically the mechanism of supra- 
transmissions observed in a boundary driven discrete nonlinear Schrodinger equation 
describing electromagnetic fields in waveguide arrays. We have shown that the source 
of the phenomenon is the presence of a saddle-node bifurcation between a stable 
discrete soliton and an unstable one. We have shown as well numerically that the 
unstable one can exhibit a different dynamics, sensitive to the perturbation. We 
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therefore argue that it might be possible to propose it as a weak signal light detector. 
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